% function Print_Equ_NK

%% Declare Variables, Parameters, and Steady States
% Variables
syms        Eta_Z Eta_Zp Eta_M Eta_Mp ...
            N w mc Y Yp Profit PAC L C ...
            Pir Pirp  ...
            T TAU B Lag_B Lag_Bp ir Lag_ir Lag_irp rr;
varlist=   [Eta_Z Eta_Zp Eta_M Eta_Mp ...
            N w mc Y Yp Profit PAC L C ...
            Pir Pirp ...
            T TAU B Lag_B Lag_Bp ir Lag_ir Lag_irp rr];
% Parameters
syms        PP_RHO_Z PP_RHO_M ...
            PP_ALPHA PP_EPSILON PP_THETA ...
            PP_Taylor_Pi PP_Taylor_ir;
pplist=    [PP_RHO_Z PP_RHO_M ...
            PP_ALPHA PP_EPSILON PP_THETA ...
            PP_Taylor_Pi PP_Taylor_ir];
% Steady States
syms        SS_Z SS_N SS_w SS_mc SS_Y SS_Profit SS_PAC ...
            SS_C SS_L ...
            SS_Pir SS_ir ...
            SS_TAU SS_B SS_T SS_G ;
sslist=    [SS_Z SS_N SS_w SS_mc SS_Y SS_Profit SS_PAC ...
            SS_C SS_L ...
            SS_Pir SS_ir ...
            SS_TAU SS_B SS_T SS_G];
        

%% Equations
%--------------------------------------------------------------------------
% Note: the aggregate quantities Y,C,G,T,Bh,Bg,PAC,Profit
%       are all normalized by SS_w*SS_N. w and N are in their original
%       levels.
%--------------------------------------------------------------------------
e       =   sym('e%d',[15,1]);
%--------------------------------------------------------------------------
% Dynamics of Exogenous Variable (2)
%--------------------------------------------------------------------------
id      =   1;
e(id)   =   -Eta_Zp+PP_RHO_Z*Eta_Z;
id      =   id+1;
e(id)   =   -Eta_Mp+PP_RHO_M*Eta_M;
%--------------------------------------------------------------------------
% Firms' Decisions (5)
%--------------------------------------------------------------------------
% Production Function 
id      =   id+1;
e(id)   =   -SS_Y*exp(Y) ...
            + SS_Z*exp(Eta_Z)*( SS_N*exp(N) )^(1-PP_ALPHA);
% Marginal Cost 
id      =   id+1;
e(id)   =   -SS_mc*exp(mc) ...
            +(SS_N*SS_w)*exp(N+w)/(SS_Y*exp(Y))/(1-PP_ALPHA);  
% Price Setting (Phillips Curves) 
id      =   id+1;
e(id)   =   -(Pir ) ...
            +1/(1+SS_Pir)*PP_EPSILON/PP_THETA*( mc ) ...
            +(1+SS_Pir)/(1+SS_ir)*(Pirp);
% Price Adjustment Costs 
id      =   id+1;
e(id)   =   -( PAC ) ...
            + 0;
% Profit 
id      =   id+1;
e(id)   =   -( SS_Profit+Profit ) ...
            + SS_Y*exp(Y)-exp(N+w)*SS_N*SS_w- ( SS_PAC+PAC );

%--------------------------------------------------------------------------
% Government (4)
%-------------------------------------------------------------------------- 
% Monetary Policy (1)
id      =   id+1;
e(id)   =   -ir ...
            +PP_Taylor_ir*Lag_ir ...
            +PP_Taylor_Pi*Pir ...
            +Eta_M;
% Fiscal Policy (3)
id      =   id+1;
e(id)   =   -( SS_T+T + (SS_B+Lag_B)/(1+SS_Pir+Pir) + SS_G ) ...
            +(SS_TAU+TAU)*( SS_w*SS_L*exp(w+L) ) ...
            +(SS_B+B)/(1+SS_ir+ir);
id      =   id+1;
% e(id)   =   -(SS_TAU+TAU)+SS_TAU;
e(id)   =   -(SS_B+B)+SS_B;
id      =   id+1;
e(id)   =   -(SS_T+T)+SS_T;

%--------------------------------------------------------------------------
% Financial Market Returns (1)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -rr+ir-Pirp;

%--------------------------------------------------------------------------
% Market Clearing Conditions (1)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -SS_N*exp(N)+SS_L*exp(L);

%--------------------------------------------------------------------------
% Lagged Variables (2)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -Lag_Bp+B;
id      =   id+1;
e(id)   =   -Lag_irp+ir;
%% Print the File
fid  =  fopen('Equ_NK.m','w+');
% Print the function head
fprintf(fid,'function [varargout]=Equ_NK(PP,SS,order)\n\n');

fprintf(fid,'%%');
for i=1:length(varlist)
    if mod(i,8)==0
        fprintf(fid,['''',char(varlist(i)),'''',',...\n %%']);
    else
        fprintf(fid,['''',char(varlist(i)),'''',',']);
    end
end

fprintf(fid,'\n \n ');
% Print the Structural Parameters
for i=1:length(pplist)
    sp  =   char(pplist(i));
    fprintf(fid,[sp,'=',strrep(sp,'PP_','PP.'),';\n']);
end
fprintf(fid,'\n ');
% Print the Steady State Values
for i=1:length(sslist)
    ss  =   char(sslist(i));
    fprintf(fid,[ss,'=',strrep(ss,'SS_','SS.'),';\n']);
end
fprintf(fid,'\n ');
% Print the Input Deviation Values
for i=1:length(varlist)
    fprintf(fid,[char(varlist(i)),'=0;\n']);
end
% Print the 0-order Equation
fprintf(fid,'if order==0\n');
fprintf(fid,['res=ones(',num2str(length(e)),',1);\n']);
for i=1:length(e)
    fprintf(fid,['res(',num2str(i),')=',char(e(i)),';\n']);
end
fprintf(fid,'varargout{1}=res;\n');
% Print the 1-order Equations
fprintf(fid,'elseif order==1\n');
for i=1:length(varlist)
    fprintf(fid,['fir_ord.',char(varlist(i)),'=[']);
    for j=1:length(e)
        der     =   jacobian(e(j),varlist(i));
        fprintf(fid,char(der));
        if j<length(e)
            fprintf(fid,';');
        end
    end
    fprintf(fid,'];\n');
end
fprintf(fid,'varargout{1}=fir_ord;\n');

fprintf(fid,'\n end\n');

fclose(fid);